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The  fundamental  processes  in  shock-induced  instabilities  of  materials  remain  obscure,  particularly  for 
detonation  of  energetic  materials.  We  simulated  these  processes  at  the  atomic  scale  on  a  realistic  model  of 
a  polymer-bonded  explosive  (3,695,375  atoms/cell)  and  observed  that  a  hot  spot  forms  at  the  nonuniform 
interface,  arising  from  shear  relaxation  that  results  in  shear  along  the  interface  that  leads  to  a  large  temperature 
increase  that  persists  long  after  the  shock  front  has  passed  the  interface.  For  energetic  materials  this  temperature 
increase  is  coupled  to  chemical  reactions  that  lead  to  detonation.  We  show  that  decreasing  the  density  of  the 
binder  eliminates  the  hot  spot. 
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The  interaction  of  shock  waves  with  nonuniform  interfaces 
plays  an  essential  role  in  the  interfacial  instabilities  in  iner¬ 
tial  confinement  fusion  (ICF),  in  shock-induced  Richtmyer- 
Meshkov  instabilities  (RMIs),  and  in  detonation  in  hetero¬ 
geneous  polymer-bonded  explosives  (PBXs).  For  detonation, 
it  is  generally  accepted  that  hot  spots  form  during  the 
development  of  instabilities  as  shock  waves  pass  through  the 
interface  or  other  defects.1-4  Despite  numerous  experimental 
and  theoretical  studies,  the  fundamental  processes  involved 
remain  controversial.  This  is  due  to  the  complex  environment 
and  coupling  of  thermal,  chemical,  and  mechanical  degree  of 
freedom,  which  is  extremely  difficult  to  unravel  experimen¬ 
tally.  It  has  also  been  very  difficult  theoretically  to  include  the 
reactive  processes  involved  and  yet  cover  the  enormous  size 
and  time  scales  intrinsic  to  the  phenomena. 

To  discover  the  origin  of  shock-induced  hot-spot  formation, 
we  carried  out  reactive  dynamics  (RD)  using  the  ReaxFF 
reactive  force  field  on  a  realistic  model  of  a  real  polymer- 
bonded  explosive  PBX  N-106.  Energetic  materials  (EMs) 
are  essential  for  applications  ranging  from  rocket  engines,  to 
building  and  dam  construction,  and  to  armaments.  Generally 
the  EM  is  bound  together  in  a  matrix  of  polymer  elastomers  to 
form  the  PBX  that  can  be  molded  into  various  shapes,  while 
providing  some  control  in  resisting  unintentional  detonation 
due  to  shocks  or  friction.  Unfortunately,  current  generations 
of  PBXs  are  sensitive  to  accidental  detonation,  despite  many 
attempts  to  control  the  safety  by  improvements  in  materials 
and  manufacturing  practices.  Here  we  use  the  ReaxFF  reactive 
force  field  to  examine  the  effect  of  shocks  on  realistic  models 
of  polymer-bonded  explosives,  where  we  use  these  simulations 
to  extract  the  mechanism  of  hot-spot  formation.  Then,  based  on 
this  model,  we  predict  how  to  change  the  system  to  reduce  the 
hot  spot,  and  carry  out  simulations  to  validate  this  prediction. 

ReaxFF  has  now  been  established  to  provide  nearly  the 
accuracy  of  quantum  calculation  in  the  various  reaction 
barriers  and  rates  while  providing  a  computational  efficiency 
nearly  that  of  ordinary  molecular  dynamics  (MD)  with 
ordinary  force  fields  (FFs),  enabling  us  to  study  the  complex 
processes  involved  in  interfacial  instabilities  at  the  atomic 
scale,  providing  insights  on  such  phenomenon.  Thus  ReaxFF 
provides  the  possibility  of  realistic  simulations  to  probe  the 


atomistic  mechanism  controlling  detonation,5-12  and  leads 
to  an  accurate  description  of  the  complex  chemistry  of 
cyclotrimethylene  trinitramine  (RDX)  under  shock-loading 
conditions,6  and  similar  calculations  on  polyethylene  (PE)  and 
poly(4-methyl-l-pentene)  polymer  lead  to  good  agreement 
with  the  experimental  Hugoniot.9 

It  is  generally  accepted  that  detonation  of  PBXs  is  initiated 
at  hot  spots,  but  the  mechanism  responsible  for  hot-spot 
ignition  is  not  clear.13  Previous  simulations14-17  in  relatively 
small  molecular  crystals  and  atomic  crystals  containing  voids 
postulated  that  void  collapse  and  nanojets  in  voids  are  the 
source  of  hot-spot  formation. 

In  order  to  develop  an  atomic-scale  understanding  of  the 
nature  of  hot-spot  formation  of  PBX  materials  we  used  ReaxFF 
RD  simulations  to  examine  the  effect  of  shocks  on  the  materials 
configuration  in  Fig.  1,  which  involves  a  sawtooth  interface 
between  the  RDX  and  the  polymer  [hydroxyl-terminated 
poly  butadiene  (HTPB)].  This  system  is  54  nm  thick  in  the 
shock  direction  (x)  with  a  period  of  27  nm  along  the  sawtooth 
direction  (y),  and  is  uniform  with  a  periodic  length  of  25  nm 
in  the  z  direction,  leading  to  a  simulation  cell  with  3695375 
independent  atoms.  For  shock  velocities  of  2.5  and  3.5  km/s 
it  takes  ~10  ps  for  the  shock  wave  to  traverse  the  interface. 
Such  a  simulation  cell  and  time  interval  is  far  beyond  the 
capabilities  of  current  quantum  mechanics  (QM)  methods 
but  can  be  treated  using  the  QM-based  ReaxFF  reactive 
force  field. 

Our  simulations  model  the  PBXN-106  material,  which 
consists  of  cyclotrimethylene  trinitramine  (RDX)  crystals 
bound  together  using  HTPB-  and  isophorone  diisocyanate 
(IPDI)-based  polyurethane  rubber.  HTPB  (see  Fig.  SI  of 
the  Supplemental  Material18)  is  a  copolymer,  containing 
80%  trans-  and  cis-butadiene  monomers  with  20%  1,2-vinyl 
monomers,  with  the  chains  terminated  by  hydroxyl  groups. 
Our  MD  simulations  use  100%  trans  monomers,  leading  to 
a  molecular  weight  of  1518  Da  per  chain.  A  plasticizer  agent 
di-octyl  adipate  (DOA)  is  added  (40%  by  weight)  to  the  binder 
to  improve  casting  properties.  To  determine  the  conformational 
state  of  the  polymer  at  300  K  and  1  atm  pressure,  we  used 
the  cohesive  energy  density  (CED)  Monte  Carlo  method19,20 
designed  to  predict  accurate  cohesive  energy  densities  of 
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FIG.  1.  (Color  online)  Snapshot  of  PBX  during  shock  loading  at 
Up  =  2.5  km/s  (for  6.0  ps).  The  shading  is  based  on  the  total  slip  in 
angstroms.  This  system  is  54  nm  thick  in  the  shock  direction  (x),  with 
a  period  of  27  nm  along  the  sawtooth  direction  (y),  and  is  uniform 
with  a  periodic  length  of  25  nm  in  the  z  direction.  The  arrow  indicates 
shock  direction. 


polymers  (CED  combines  repeated  cycles  of  temperature 
annealing  and  quenching  simultaneous  with  density  annealing 
and  quenching  to  ensure  that  polymer  chains  are  fully 
equilibrated).  We  used  CED  to  build  ten  independent  samples 
of  each  polymer  within  a  unit  cell  containing  two  independent 
molecular  chains  and  eight  DOA  molecules.  We  chose  the 
mass  ratio  of  binder  over  DOA  to  be  1.03,  which  leads 
to  the  optimum  viscosity  and  compression  strength.21  Each 
chain  contains  ten  HTPB  repeat  units  connected  via  one  IPDI 
crosslinking  molecule  to  four  terminal  HTPB  repeat  units  per 
end.  The  initial  conformation  was  based  on  rotational  isomeric 
states  (RISs)  statistics  to  determine  the  dihedral  angle  values 
for  each  chain  and  its  repeating  units.  The  torsional  barriers 
for  the  RIS  and  the  charges  (Mulliken)  were  determined  [from 
Becke  three-parameter  Lee- Yang-Parr  (B3LYP)  hybrid  density 
functional  calculations  using  the  6-311  G**  basis  set]  on 
a  polybutadiene  trimer.  The  same  level  of  QM  calculations 
was  also  performed  on  the  IPDI  and  DOA  to  determine  the 
charges  and  structures  for  the  binder  model.  These  QM  results 
and  model  preparation  procedure  are  provided  as  part  of  the 
Supplemental  Material.18 

To  simulate  the  propagation  of  shock  waves  in  the  PBX, 
we  impacted  the  two-dimensionally  (2D)  periodic  PBX  slab 
onto  a  Lennard  Jones  9-3  wall.  The  shock  direction  (x) 
was  set  to  be  finite.  The  slab  was  thermally  equilibrated  at 
T  =  300  K  for  4  ps  and  then  propelled  at  the  desired  impact 
velocity  by  adding  a  corresponding  translational  velocity  on 
top  of  the  thermal  component.  This  procedure  produces  a  shock 
wave  with  an  average  particle  velocity  (Up)  equal  to  the  impact 
velocity,  propagating  into  the  slab  from  the  RDX  side.  Here 
we  explore  the  shock  response  of  the  nonuniform  interface  in 
the  PBX  at  Up  =  2.5  and  3.5  km/s  via  adiabatic  MD  to  follow 
the  dynamic  evolution  of  the  system,  using  a  time  step  of  0. 1  fs 
to  integrate  the  equations  of  motion,  for  up  to  10  ps. 

Considering  the  heterogeneities  normal  to  the  shock  direc¬ 
tion  (x),  we  partition  the  simulation  cell  into  fine  bins  along 
the  x,y  directions  and  average  all  physical  properties  of  interest 
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within  each  bin,  including  stress  (Oij)  and  temperature  ( T ) 
profiles  in  the  xy  plane  at  different  stages  of  compression 
(i,j  =1,2,  and  3,  corresponding  to  x,  y,  and  z,  respectively.)  The 
atomic  stress  is  obtained  from  the  atomic  virial.22  The  stress 
calculation  is  based  on  the  center  of  mass  (COM)  of  each 
molecule  to  eliminate  binning  effects.  The  bin  width  is  ~1  nm 
x  1  nm  in  the  xy  plane.  The  von  Mises  shear  stress  r  is  also 
computed  as  the  maximum  shear  stress,  as  2r  =  a xx  —  (ayy  + 
cfzz)/ 2,  under  the  assumption  of  ayy  =  o  zz,  which  is  true  in  our 
simulations. 

To  track  the  chemical  processes  as  the  shock  wave  propa¬ 
gates  through  the  PBXs,  we  analyze  the  fragments  resulting 
from  the  corresponding  ReaxFF  MD  trajectories  based  on 
a  time- windowed  average  of  bond  orders.  Combining  the 
fragment  analysis  results  and  the  two-dimensional  binning 
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FIG.  2.  (Color  online)  Time  evolution  of  temperature  and  von 
Mises  shear  stress  for  Up  =  2.5  km/s  at  various  stages  during  the 
shock  compression.  The  solid  line  represents  the  interface  position 
between  the  RDX  and  the  binder.  Shock  wave  propagates  from  left 
(RDX)  to  right  (binder),  (e)  4.5  ps,  shock  wave  has  not  propagated 
to  polymer;  (a),(f)  5.0  ps,  shock  wave  just  reaches  the  interface; 
(b),(g)  5.5  ps,  adiabatic  shear  localization  in  a  small  triangle 
interface  region;  (c),(h)  6.0  ps,  shock  wave  has  passed  through 
the  hot-spot  region;  (d)  8.0  ps,  shock  wave  has  passed  through  the 
interface. 
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analyses,  we  examined  the  chemistry  in  the  xy  plane  at 
different  stages  of  shock  loading.  To  examine  the  atomic  local 
shear,  we  use  maximum  relative  displacement  (MRD), 23,24 
which  is  defined  as  s ,•  =  x^-  —  X;7-:lx/7-  —  Xf-7lmax.  Here 
the  Xij  and  X/7  vectors  correspond  to  the  difference  between 
the  current  and  reference  configurations  for  atoms  i  and  j 
(slipped  neighbors  of  atom  i),  respectively.  The  reference 
configurations  are  taken  to  be  the  preshock  structures. 

Figures  2(a)-2(d)  show  the  2D  profile  of  temperature  T 
and  shear  stress  r  at  various  stages  of  compression  for  Up  = 
2.5  km/s.  The  temperature  at  the  angular  tip  interface  region 
increases  by  ~400  K  compared  to  other  shocked  regions,  once 
the  shock  wave  has  passed  through,  leading  to  formation  of 
a  hot  spot.  The  hot-spot  region  expands  into  the  explosive 
along  the  interface  as  the  shock  wave  continues  to  propagate. 
The  shock  wave  converges  as  it  propagates  in  the  binder  due 
to  the  wave  refraction,  as  it  passes  through  the  nonuniform 
interface.  This  shock-wave  convergence  leads  to  a  small, 
second  highest  temperature  region  ~50  K  higher  (see  Fig.  S5 
in  the  Supplemental  Material18)  to  the  right  of  the  hot-spot 
region  in  the  binder.  This  hot  spot  also  appeared  at  higher 
shock  velocities.18 

To  understand  the  mechanism  underlying  the  hot-spot 
formation  found  in  our  study,  we  traced  the  von  Mises  shear 
stress  at  different  stages.  Upon  shock  arrival  to  the  interface, 
we  find  reduction  of  shear  stresses  at  the  angular  tip  regions 
as  manifested  in  Fig.  2(g).  This  arises  from  rapid  relaxation 
of  the  stress  at  the  tip  by  the  more  compressible  polymer, 
possibly  combined  with  the  convergence  of  the  shock  wave 
at  the  tip,  interfacial  friction  between  RDX  and  binder,  and 
perhaps  the  collapse  of  tiny  nanovoids  at  the  angular  tip.  This 
shear  localization  at  the  angular  tip  leads  to  nucleation  of  the 
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hot  spot.  The  stress  decreases  rapidly  in  the  hot-spot  region, 
leading  to  large  local  shear  deformations  as  shown  in  Fig.  2(h). 
This  small  hot-spot  region  expands  rapidly  as  the  shock 
wave  propagates,  causing  additional  shear  localization  near 
the  interface  in  the  explosive  and  an  additional  increase  in  the 
local  temperature.  Since  the  time  scale  is  only  picoseconds,  we 
consider  this  to  be  an  adiabatic  shear  localization  process.  The 
expansion  in  the  shear  localized  region  is  depicted  in  Figs.  2(g) 
and  2(h).  The  MRD  analysis  in  Fig.  1  shows  that  the  regions 
with  maximum  local  deformation  correspond  to  the  maximum 
local  shear  stress  relaxation,  and  to  the  hot  spot.  Since  the 
shock  wave  propagates  from  high-density  (RDX)  explosives  to 
a  low-density  (binder)  polymeric  matrix,  the  reflective  waves 
from  the  interface  will  decrease  the  stress  in  the  explosive, 
which  tends  to  suppress  the  hot-spot  formation  as  discussed 
below.  Thus  we  find  that  the  hot  spot  forms  in  our  system  as  a 
compromise  between  shear  localization  and  reflective  wave  re¬ 
laxation  in  the  angular  tip  and  interface  regions.  Our  atomically 
based  hot-spot  ignition  mechanism  provides  an  understanding 
that  differs  from  previously  proposed  mechanisms25  involving 
adiabatic  compression  of  a  trapped  gas,  cavity  collapse 
including  viscous  or  plastic  heating  of  the  surrounding  matrix 
material,  friction  between  the  impacting  surfaces,  among  other 
phenomena. 

In  order  to  follow  the  chemical  processes  as  the  shock  wave 
propagates  through  the  interface  of  the  RDX  and  the  binder, 
we  plot  the  distribution  of  NO2  production  in  the  xy  plane 
since  NO2  dissociation  is  the  main  chemical  product  for  Up 
below  4  km/s  (see  Ref.  6).  Figure  3  shows  the  time  evolution 
of  NO2  for  Up  =  2.5  and  3.5  km/s.  We  find  few  reactions 
for  Up  =  2.5  km/s,  except  those  occurring  in  the  hot-spot 
region.  The  plastic  shear  deformation  in  the  hot-spot  region 


(a)  Up  =  2.5  km/s,  time  =  6  ps  (b)  Up  =  2.5  km/s,  time  =  8  ps 


(c)  Up  =  3.5  km/s,  time  =  4.5  ps  (d)  Up  =  3.5  km/s,  time  =  6.5  ps 


FIG.  3.  (Color  online)  2D  distribution  and  time  evolution  of  chemistry  (N02  formation)  for  various  Up.  (a)  U p  =  2.5  km/s  at  t  =  6.0  ps: 
The  chemical  reaction  has  occurred  mainly  in  the  hot-spot  region,  (b)  Up  =  2.5  km/ s  at  t  =  8.0  ps:  Additional  chemical  reactions  have  occurred 
and  a  small  amount  of  N02  has  diffused  to  the  polymer  region  (the  right  region  of  the  solid  line  in  the  xy  plane),  (c)  Up  =  3.5  km/s  at  t  =  4.5 
ps:  The  chemical  reaction  occurs  in  the  shocked  regions  but  much  more  occur  in  the  hot-spot  region,  (d)  Up  =  3.5  km/s  at  t  =  6.5  ps:  We  see 
diffusion  of  N02  into  the  polymer  region. 
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FIG.  4.  (Color  online)  The  dependence  of  hot-spot  formation  (temperature  and  shear  stress)  on  the  density  of  the  polymer  binder  at  a  shock 
velocity  of  3.5  km/s.  The  middle  (0.95  g/cm3)  is  normal  density.  Thus  using  of  a  polymer  with  half  the  normal  density  prevents  hot-spot 
formation.  This  is  due  to  the  improved  interface  impedance  which  allows  reflective  rarefaction  wave  relaxation.  In  contrast,  increasing  the 
binder  density  enhances  the  hot  spot  due  to  the  reflective  shock  waves  reflected  back  into  the  RDX. 


produces  enough  energy  for  N-N  bond  excitations,  leading  to 
increased  NO2  formation.  The  same  phenomenon  occurs  for 
JJp  =  3.5  km/s.  The  chemical  reactions  occurring  in  the  hot¬ 
spot  region  release  additional  energy,  expanding  it  to  induce 
even  more  shear  deformation.  Under  the  coupling  effects  of 
local  shear  relaxation,  increased  chemical  reactions,  increased 
local  temperature,  and  the  reflective  wave  relaxation,  the  hot 
spot  expands  until  detonation  (or  at  lower  input  shock  velocity 
it  can  disappear). 

Our  hot-spot  ignition  mechanism  suggests  that  use  of  a 
binder  with  lower  density  would  lead  to  decreased  hot-spot 
formation  and  hence  a  less  sensitive  EM  for  propulsion  and 
explosives  applications.  To  test  this  idea  we  changed  the 
binder  density  to  half  the  original  value  (0.95  g/cm3),  to 
approximately  twice  the  value  by  scaling  the  atomic  mass.  We 
carried  out  shock  simulations  at  Up  =  3.5  km/s.  To  analyze 
the  hot-spot  behavior  we  display  a  2-nm  segment  in  the  y 
direction  for  a  one-dimensional  binning  analysis  during  shock 
compression.  Figure  4  shows  the  temperature-time  and  shear 
stress-time  diagrams  from  3.6  ps,  at  which  the  shock  wave 
has  not  reached  the  interface  region.  The  temperature-time 
diagram  shows  that  the  hot  spot  is  dramatically  increased 
for  the  high-density  binder,  but  it  disappears  for  the  half¬ 
density  case  (0.48  g/cm3).  Indeed,  this  low-density  binder 


suppressed  hot-spot  formation  due  to  the  relaxation  of  the  hot 
spot. 

Previous  studies26,27  showed  that  a  second  shock  wave 
(double  shock)  or  a  rarefaction  wave  will  be  reflected  back 
into  the  original  medium  as  the  shock  wave  travels  through 
interface.  These  refraction  wave  patterns  are  determined  by 
the  shock  impedance,  which  is  a  function  of  the  incident  wave 
angle,  the  sound  speed  at  the  shocked  materials  interfaces, 
and  the  stresses  in  the  two  regions.  As  shown  in  the  normal 
stress-time  behavior  (crn-t,  see  Fig.  S6  of  the  Supplemental 
Material18),  a  high-density  polymer  matrix  leads  to  a  reflected 
shock  wave  due  to  the  higher  shock  impedance  at  the  interface, 
while  the  low-density  polymer  matrix  leads  to  a  rarefaction 
wave.  The  von  Mises  shear  stress  diagram  in  Fig.  4  shows  the 
same  character:  A  high  shear  stressed  reflective  wave  forms 
for  the  high-density  case,  a  modest  one  for  normal  density, 
and  complete  disappearance  for  the  low-density  case  (here 
the  reflective  wave  relaxes  the  system,  avoiding  the  hot  spot). 
These  results  indicate  that  the  hot  spot  arises  from  the  complex 
variation  of  shock  impedance  at  this  nonuniform  interface.  We 
also  note  from  the  x-t  diagram  that  the  shock  wave  propagates 
at  a  slower  speed  through  the  high-density  binder  matrix. 

These  results  suggest  that  a  critical  design  strategy  for 
developing  insensitive  EM  for  propulsion  and  explosives  is 
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to  use  a  binder  with  significantly  lower  mechanical  impedance 
than  that  of  EM.  Thus,  merely  reducing  the  density  of  the 
polymer  by  a  factor  of  2  eliminates  the  hot  spot. 

Summarizing,  we  used  reactive  dynamics  (RD)  to  examine 
the  mechanical,  chemical,  and  thermal  response  of  mechan¬ 
ically  shocked  polymer-boned  explosives  (PBXs)  by  using 
a  realistic  model  of  a  nonuniform  interface  with  3695375 
independent  atoms.  We  observe  that  a  hot  spot  develops 
at  the  first  point  at  which  the  shock  wave  encounters  the 
lower-density  polymer  at  the  nonuniform  interface  between  the 
explosive  and  elastomer.  Our  analysis  shows  that  the  hot  spot 
arises  from  adiabatic  shear  localization,  which  leads  to  a  large 
temperature  increase  that  persists  long  after  the  shock  front 
has  passed  the  interface,  followed  for  energetic  materials,  by 
coupling  to  chemical  reactions  (NO2  and  HONO  dissociation) 
that  ultimately  leads  to  detonation. 

We  should  emphasize  the  sequence  of  events  here: 
(a)  Shock  at  the  asperity  of  a  nonuniform  interface  of 
heterogeneous  materials  leads  to  shear  relaxation  that  causes 
shear  stresses  along  the  nonuniform  interface,  (b)  These  shear 
stress  relaxations  result  in  energy  deposition  at  the  interface 
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that  leads  to  a  local  significant  increase  in  temperature,  (c)  This 
increase  in  temperature  leads  to  bond  breaking  in  energetic 
materials,  which  leads  eventually  to  detonation  (which  may  be 
premature). 

For  nonenergetic  materials  we  also  expect  the  first  two 
phenomena  of  local  shear  stress  and  temperature  increase. 
Thus  for  nonenergetic  materials  we  can  expect  that  the  hot 
spot  may  cause  local  melting  and  defect  generation  that  would 
accelerate  failure  and  perhaps  instabilities  (e.g.,  Richtmyer- 
Meshkov). 

These  results  suggest  that  a  critical  design  strategy  in  devel¬ 
oping  insensitive  EM  for  propulsion  and  explosives  is  to  use  a 
binder  with  a  much  lower  density  (1/3  that  of  the  explosive). 
We  tested  this  by  modifying  the  polymer  system  to  have  half 
the  normal  density  (1/3  that  of  the  explosive)  and  showed  that 
no  hot  spot  develops  using  the  same  shock  conditions. 

All  computations  were  carried  out  on  the  Army  HPC 
systems;  we  thank  Betsy  Rice  and  Larry  Davis  for  assistance. 
Personnel  were  supported  by  ARO  (W911NF-05-1-0345  and 
W91  INF-08-1-0124)  and  by  ONR  (N000 14-09- 1-0634). 
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